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Abstract. We obtain approximations for the CDM particle trajectories starting from La- 
grangian Perturbation Theory. These estimates for the CDM trajectories result in approxi- 
mations for the density in real and redshift space, as well as for the momentum density that 
are better than what standard Eulerian and Lagrangian perturbation theory give. For the 
real space density, we find that our proposed approximation gives a good cross-correlation 
(> 95%) with the non-linear density down to scales almost twice smaller than the non-linear 
scale, and six times smaller than the corresponding scale obtained using linear theory. This 
allows for a speed-up of an order of magnitude or more in the scanning of the cosmologi- 
cal parameter space with N-body simulations for the scales relevant for the baryon acoustic 
oscillations. Possible future applications of our method include baryon acoustic peak recon- 
struction, building mock galaxy catalogs, momentum field reconstruction. 



1 Introduction 

Understanding the evolution of the large-scale structure (LSS) of the universe is a crucial 
ingredient of present-day cosmology. Because of the larger accessible volume, the importance 
of the LSS will only grow in the future as it has the potential to become an even more accurate 
cosmological probe than the Cosmic Microwave Background (CMB). Thus observations of the 
LSS alone or combined with CMB observations could provide even stronger constraints on the 
equation of state of dark energy, modifications to gravity and primordial non-Gaussianities. 

On the theoretical side, the largest scales in the universe (> 100 Mpc at redshift z = 0) 
can be treated perturbatively (e.g. [1]). However, the perturbative treatment of the LSS 
breaks down when the overdensity (the standard parameter controlling the perturbative 
expansions) becomes of order 1. This corresponds to the mildly non-linear (MNL) regime: 
scales < 50 Mpc at z = 0. The dynamics at these scales is important as it introduces 
percent-level shifts to the acoustic peak caused by the Baryon Acoustic Oscillations (BAO) 
in the matter correlation function [2]. On the observational side, such an accuracy in the 
determination of the position of the acoustic peak is required by ongoing and future surveys 
of LSS, such as BOSS 1 and WFIRST 2 . Such surveys target biased tracers of the underlying 
density field in redshift space. Thus, to make contact with observations, what is required 
is a good understanding not only of the MNL regime, but of the non-linear (NL) regime as 
well, since finger-of-god (FoG) effects from NL structures, such as galaxy clusters, leak into 
the MNL scales. 

Most theoretical studies of structure formation in the MNL and NL regimes are currently 
done through numerical simulations. However, computational time is still a very important 
constraint on simulations. Sampling variance requires one to perform numerous simulations 
of different realizations of the same cosmology; or alternatively simulate large volumes. This 
prevents an efficient sampling of the cosmological parameter space (e.g. [3]), and thus, it 
is still difficult to get a complete handle on the uncertainties that would arise from mea- 
surements at the MNL scales relevant to the BAO, and the NL scales relevant to e.g. weak 
lensing and galaxy clustering. 

Having an accurate model for the (M)NL regimes of structure formation is therefore 
important for the utilization of observations of galaxy clustering. One important application 
of such a model will be for improving the reconstruction techniques [4, 5] of the baryon 
acoustic peak from observations. One can argue that reconstruction pulls the information 
stored in the higher-re-point functions and puts it back into the power spectrum, and thus 
has been proven to reduce the errors in the determination of the position of the acoustic peak 
[4, 5] - a necessary improvement for utilizing current and future BAO observations. 

In a recent paper [6] (TZ from now on), we offered a model for the matter distribution in 
the MNL regime. We showed that through an expansion around the Zel'dovich approximation 
(ZA) [7] (or higher-order Lagrangian Perturbation Theory (LPT) [8]) one can construct 
robust approximations to the MNL density field in real space. Therefore, one can write the 
true overdensity 5 in terms of the cheap-to-compute density field in the ZA, S z , as (in Fourier 
space) : 

S = R s 5 z + 5 M c, (1.1) 



1 http://www.sdss3.org/surveys/boss.php 
2 http://wfirst. gsfc.nasa.gov/ 



where R$ is a non-random scale- and time-dependent transfer function; and 5mc (MC stand- 
ing for mode coupling) is a small random residual at MNL scales. We showed that the re- 
sulting matter power spectrum estimator reduces sample variance by more than an order 
of magnitude in the MNL regime. We found that 5 and 5 Z are strongly correlated in that 
regime, which implies a small 8mc an d an R$ which has small sample variance: a result 
which is key to the success of our treatment in TZ. 

In TZ we also showed that other methods (e.g. [9, 10]) based on expansions similar to 
eq. (1.1) but done around the linear density field, 5l, do not offer similar improvements to 
cosmic variance. We found that this is caused by the fact that 5l and 5 become decorrelated 
for scales smaller than the rms particle displacements (corresponding to a wavevector k ~ 
0.12/i/Mpc). We refer the reader to TZ for a detailed discussion of this result. 

In general, in this paper we refer to quantities Q+ which are highly-correlated to the 
corresponding non-linear quantities Q, as approximations of Q. This is because, once we 
know the transfer function Rq relating Q and Q+, we can construct approximations to Q 
by writing RqQ+. Obviously, to build unbiased estimators, one needs to include the effect 
of the MC terms: Q = RqQ+ + QmCi as done in eq. (1.1) and in TZ. However, since by 
construction Q+ and Q are strongly correlated, we have that Qmc is small and Rq has small 
sample variance and can be obtained by using only few modes. In this terminology, 5 Z is a 
good approximation of 8, while 8l is not. 

One problem with the TZ treatment is that it is not suitable for including the effects of 
Redshift Space Distortions (RSD). It is not particle-based, i.e. it constructs an approximation 
of the nonlinear density by treating it as a field, and not as generated by discrete particles 
with known velocities. Furthermore, this prevents one from extending the TZ method for 
constructing mock galaxy catalogs based on estimated particle positions and velocities. An- 
other way of saying this is that working only with the density field misses the complicated 
phase-space structure of Cold Dark Matter (CDM), which can be recovered in the particle 
picture. 

In this paper, we address this issue by concentrating on the CDM particle trajectories in 
LPT, and relating those to the trajectories obtained from N-body simulations. The reason we 
choose to work around the LPT results is again related to the discussion in TZ of the effects 
of the bulk flows, which cancel explicitly at each order in LPT. We find that using approxi- 
mations of the particle velocities and positions based on LPT results in good approximations 
for the density in real and redshift space, which allows us to go well into the MNL regime 
down to scales comparable and even smaller than the non-linear scale (Jznl ~ 0.25/i/Mpc for 
redshift of z = 0). 

In a nutshell, the method used in this paper for building an approximation of a quantity, 
Q[x,v], which is a functional of the particle positions (x) and velocities (v) is as follows. 3 
It consists of first calculating the displacement (related to the positions xlpt) and velocity 
( v lpt) fields in LPT. Those are then transformed using transfer functions obtained after 
a calibration with simulations. These transfer functions, R and R", relate the LPT dis- 
placement and velocity fields to their respective fully non- linear counterparts. The effect 
of these transfer functions is twofold: 1) To capture higher-order LPT effects, which are 
not included under a truncation of the LPT expansion; 2) To capture some of the effects 
of stream crossing, which are missed by LPT (see below). Therefore, schematically we can 
write x = x+ + XmC = R* x^pt + xmc and v = v+ + vmc = R v * vlpt + Vmc, where 



3 Here Q can stand for the matter density in real and redshift space, as well as the momentum and kinetic 
energy densities. 
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Figure 1: Shown are various matter power spectra at z = for ACDM, similar to Fig. 1 of 
[6] (errorbars are suppressed for clarity). The power spectra are divided by a smooth BBKS 
[11] power spectrum with shape parameter T = 0.15 in order to highlight the wiggles due 
to the BAO. The solid curve shows the non-linear power spectrum (i.e. the "exact" power 
spectrum obtained from N-body simulations); the part of the NL power spectrum due to 
the "memory of the initial conditions" [9] is shown with the dashed line; the power due to 
the projection of the non- linear density field on our best approximation (cf. second line of 
eq. (3.4)) is given by the dotted curve. The residual mode-coupling power is 13.5% of the 
NL power at k = 0.5/i/Mpc for z = 0. Here R$ corresponds to the transfer function relating 
the NL density and the density given by our best approximation. A vertical line denotes the 
non-linear scale. The current Hubble expansion rate in units of lOOkm/s/Mpc is given by h. 



* denotes a convolution in Lagrangian space, and quantities with subscript MC are random 
residuals. One can then construct Q+ = Q[x+,v+}. The effect of the MC terms can then 
be captured by relating Q+ to Q by using yet another transfer function (Rq) and a mode 
coupling term (Qmc)'- Q = RqQ+ + Qmc (i n Fourier space). The main benefit of first 
transforming the LPT displacement and velocity fields is that, as we will show, the phase- 
space structure of the CDM is recovered better. This means that Q+ as calculated above is 
better correlated to Q, than if Q+ was built from the uncorrected LPT quantities (such that 
Q-k = Qlpt = Q[xlpt,vlpt}) as we did in TZ. 

One may be worried that the proliferation of transfer functions compared to TZ is an 
overkill. However, those can be obtained from calibration with modest N-body runs as the 
transfer functions have rather small sample variance. Another point we would like to stress is 
that in this paper we concentrate on the CDM dynamics only, neglecting the effects of baryon 
physics. Those, however, should be possible to capture by calibrating the displacement and 
velocity transfer functions with N-body simulations, which include baryon physics. 

To demonstrate the power of the method described above, in Fig. 1 we show the NL 
power spectrum obtained from simulations, along with the power spectrum of R$5+ (without 
the MC power) obtained from our method (cf. Fig. 1 in TZ). For comparison, we include the 
power spectrum of the piece of the overdensity (Rl^l) corresponding to the memory of the 
initial conditions [9] . Thus, the residual mode-coupling power for the method described above 
is indeed rather small well into the NL regime - it is 13.5% of the NL power at k = 0.5/i/Mpc 



for z = 0. As we will show later on, the estimators for the NL power spectrum based on the 
method above (see also TZ) have a negligible variance in the linear regime (the variance is 
suppressed by a factor of ~ 10 5 for k ~ 0.02/i/Mpc at z = 0). The cross-correlation coefficient 
between our best 5+ and the NL 5 is > 0.95 up to k ~ 0.45/i/Mpc at z = 0. One should 
compare this number with the corresponding scale for the pure 2LPT-based approximation 
(the best approximation used in TZ): k ~ 0.37/i/Mpc. 

However, one is ultimately interested in tracers of the underlying density field. There- 
fore, in Fig. 2, we show the probability density for the distribution of errors in approximating 
the center-of-mass (CM) positions and velocities of Eulerian volume elements; as well as the 
distribution of errors in estimating the root-mean-square (rms) velocity, responsible for the 
FoG effect. We show the probabilities both for all space (1+5+ > 0) and for overdense regions 
(1 + 5+ > 1) only. We obtained Fig. 2 after Gaussian smoothing on a scale of 1.3 Mpc//i and 
therefore 1 + 5+ > 1 corresponds to regions containing approximately a galaxy mass or more. 
From the figure one can argue that our method should be robust to within a couple of Mpc 
in predicting tracer positions both in real and in redshift space (neglecting the FoG). The 
FoG effects can lead to somewhat larger errors. However, the BAO scale is determined by 
averaging the cross-correlation between numerous tracers, which implies that the precision 
in measuring the BAO scale should be better than the precision in the position of individual 
tracers. Comparing the curves for 1 + 5+ > and 1 + 5+ > 1 one can also argue that our 
method works much better in underdense regions - a result which we confirm later on. 

In Section 2 we write down our model for the particle trajectories, and their relation to 
the density in real and redshift space. In Section 3 we show the displacement and velocity 
transfer functions obtained from simulations, and show their relation (at lowest order in 
LPT) by building a simple model for their time dependence. Then, we plot the resulting 
density transfer and cross-correlation functions in real and redshift space (as well as the ones 
for the momentum and kinetic energy densities), and compare our results with the results in 
TZ. We summarize in Section 4. 

2 Expanding Perturbatively the CDM Equation of Motion 

2.1 Density in real and redshift space 

In this section we will relate the CDM particle trajectories (those can be perturbative or 
obtained from N-body simulations) to the density in real and redshift space. These relations 
will later allow us to discuss the resulting approximations for the density, starting from the 
approximations for the particle trajectories, which we will construct later on. 

We start off by writing the relation between the Lagrangian coordinates, q, of a particle 
and its Eulerian comoving coordinates, x: 

x(q,v) =Q + s(q,rj) , (2-1) 

where rj is conformal time; and s is a stochastic displacement vector field. As an illustration, 
in the Zel'dovich approximation [7], s is given by s z (q,rj) = D(r])s z (q,rjo), where D is the 
growth factor with D(rjo) = 1. Here d q • s z = —5l, where 5l is the fractional overdensity 
given by linear theory. 

The fully non-linear fractional CDM overdensity in real space is then given by 

5(x, n)= f d 3 q5 D (x - q - s) - 1 , (2.2) 
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Figure 2: Probability density, p, for the distribution of errors in determining the center-of- 
mass (CM) position/ velocity of a given volume element; as well as the error in determining 
the rms velocity (z = 0). We show the distributions after imposing the denoted cutoffs in 
the approximated density field, 5+. The approximate velocities, positions and density are 
calculated according to our best model (cf. second line of eq. (3.4)). CM velocities and 
positions for high-density regions are determined with worse precision than for voids, but 
still the error is within a couple of Mpc. The rms velocity responsible for the FoG effect is 
determined with worse precision but again within several Mpc. The rms linear displacement 
corresponding to z = is ~ 15Mpc/h, or about an order of magnitude larger than the errors 
depicted in the plot. The probabilities are obtained after Gaussian smoothing on a scale of 
1.3Mpc//i. 



where 5e> is the Dirac delta function. The integral above is simply a sum over all CDM 
streams of the number of particles falling in an Eulerian volume element around x. 

Similarly, we can write the redshift space coordinates, x s (q,rj), of a particle as x s = 
x + Su/'H = q + s + SufH, where sm is the line-of-sight (LoS) piece of the peculiar velocity, 
s. A dot represents a derivative with respect to conformal time, and T~L = a/a = aH, where 
a is the cosmological scale factor, and H is the Hubble parameter. Therefore, the fractional 
overdensity in redshift space, 6 S , is given by 



5 s (x s , 7]) = / d 3 q5 D (x s -q-s 



H 

n 



i 



S D (x-q- s) - 1 . 



(2.3) 



We can Fourier transform with respect to x s and expand with respect to the velocity to get 
(M0): 



6 s (k,rj) = Y^ 



1 



n=0 



-ik\ 



n! V H 



Tfr\k, V ), 



(2.4) 



where T!J n (fc, rf) is the Fourier transform with respect to x (not x s ) of the quantity: 

T^(x,ri)= fd 3 q5 D (x-q-s(q,ri)]s\l(q, V ) . (2.5) 



The result, eq. (2.4), matches that of [12] (up to a Fourier convention). 

Therefore, as in [12] we can see that at zero order in the velocity, 5 S is given by the 
real-space density 5. At first order it is related to the LoS momentum density; while at 
second order - to the LoS kinetic energy density. 

2.2 The CDM equation of motion 

To see how the displacement field can be expanded perturbatively, we have to start with the 
equation of motion for CDM particles: 

x + nx = -V4> , (2.6) 

where <j) is the gravitational potential sourced by 5 according to the Poisson equation. To 
write eq. (2.6), we restrict ourselves to modes well inside the horizon, and neglect the effect 
of baryons. Using eq. (2.1,2.2), the above equation can be rewritten as 

s"(q^(T)) = ~VV- 2 5 (2.7) 



~ / ^H)p eife ' M ex P ( ik ■ (*(<?, 17) - *(«, V)) ) 



using Jeans swindle in the form of 5n{k)k/k 2 = (i.e. a uniform density field generates 
zero acceleration). Here primes denote a derivative with respect to r defined as dr = drj/a. 
Unlike LPT, the equation above is exact, containing all information about stream crossing. 4 
Since we want to go beyond LPT, another route for solving eq. (2.7) is needed. 

The above equation can in principle be solved using the Born approximation starting 
from the Zel'dovich solution s z . This automatically recovers the particle interpretation of 
the Helmholtz hierarchy developed in [13] (see their eq. (8.11) with m = n — 1). However, 
solving this hierarchy still requires intensive N-body calculations (see [13]). 

Another way of thinking about eq. (2.7) is as an equation of motion for the vector 
field s in a non-linear theory, whose first order ("free theory") solution is the Zel'dovich 
approximation in which s is a Gaussian random field. However, as can be seen from eq. (2.7), 
the interaction vertices for s are infinitely many and of infinite order. Therefore, treating s 



4 To compare eq. (2.7) with LPT, we need to take the divergence (in Eulerian space) of that equation. We 
obtain: 



Tr 



1 + M v ; 



- -Pis - — - — V i 

D ~ D D 4r* det [1 + M (g) 



(2.8) 



= -p ~ -p j (2tt) 3 e ' fc(q ^ 6XP \ ik ' ^ Q ' ^ ~ S ^' V ^ 

where (M)ij — d qi Sj is the deformation tensor. The sum is restricted over all streams originating from q 
which arrive at x = q + s(q,rj) at time rj. LPT (e.g. [8]) is obtained by writing the density entering above in 
the single-stream approximation. That implies restricting the sum on the first line of eq. (2.8) to q — q, and 
replacing Si(q) — $i(q) on the second line with (q — q)jMji(q). 



analytically in the interaction picture and not using the single-stream approximation for the 
density is still an immensely complicated task. 

Therefore, in the rest of this paper, we will try to solve eq. (2.7) by guessing a solution 
for s, which order by order is given by the LPT result, but convolved with a transfer function 
to be obtained by a calibration with simulations. The effect of the transfer functions will be 
twofold: 1) To capture higher-order LPT effects, which are not included under a truncation 
of the LPT expansion; 2) To capture some of the effects of stream crossing, which are missed 
by LPT (see footnote 4). 

2.3 Expanding the displacement field 

In this section we will write down the approximations for the CDM trajectories following 
the prescription at the end of Section 2.2. We build a linear model for the fully non-linear 
displacements and velocities starting from the LPT results (in analogy with the linear model 
for the density we built in TZ). For simplicity, we restrict ourselves to second order LPT. 
Thus, we can write (in Fourier space with respect to q) : 

s(k, rj) = s+(k, rj) + s MC (k, rj) (2.9) 

= R z (k,T])s z (k,r]) + R2LPT(k,r])s 2 LPT(k,r]) + s M c(k,ri) 

v(k, rj) = x(k, rj) = s(k, rj) = v+(k, rj) + v MC (k, rf) (2.10) 

= Ki k i V)$z(k, rj) + R v 2LPT {k, r])s 2L PT{k, rj) + v M c(k, rj) , 

where s+ = s — Smc an d v+ = v — Vmc are the displacement and velocity approximations 5 , 
respectively; a dot denotes a derivative with respect to rj; and k above denotes wavevectors 
corresponding to q (not x); R z /2LPT(k, v) is the transfer function acting on the displacement, 
s z/2LPT(k, rj), obtained at the corresponding order in LPT. 6 We denote with smc the residual 
("mode coupling") displacement, defined such that {8Mc(k,r])s*, 2LpT (k,r])) = 0, where 
angular brackets denote ensemble averages. The displacement transfer functions are given 
by a least squares fit in Lagrangian space: 7 

(s i (k,r ] ) S * i (k,r,)) 
R z (k,rj) = (no summation) 

(s*,t(MK,i(M)) 

(si(k, rj)s2 LPT { (k, rj)) 

R2LPT{k,rj) = — - -! — (no summation) , (2.11) 

{S2LPT,i{k^) s 2LPTA k ^)) 

where we used that {s z S2lpt) = in the absence of primordial non-Gaussianities, since s z is 
first order in 5l, while S2lpt is second order. Obviously, the transfer functions above depend 
on the order at which we truncate the LPT expansion. However, since {s z S2lpt) = 0, R z 
and R v z are the same whether we truncate at first or second order in LPT. Following TZ, 
in order to make the above transfer functions unbiased, when numerically calculating the 
ensemble averages one has to average over all available modes (i.e. both over all simulation 



5 Note that we do not need extra transfer functions Rq (see eq. (2.12)) when Q stands for s or v, because 
of the linearity of our model (2.9, 2.10), which means that any Rq can be absorbed in the R's and R v, s. 

6 The lowest order in LPT is the the ZA (subscripts z), and the next is 2LPT (subscripts 2LPT). 

7 Note that the displacement and velocity transfer functions should in principle be tensors. Thus, for 
example, the term R z s Zt i should really be R z ,ijS z ,j- However, due to isotropy R z ,ij is proportional to the 
identity matrix, which allows us to replace it with the scalar transfer function, R z . One can make the same 
argument for the rest of the transfer functions. 



boxes and over all k). As long as s and s z / 2 lpt are wen correlated, Smc is small, and the 
-R's have small sample variance, and can be obtained from relatively few modes in the same 
way we obtained the density transfer function in TZ. The corresponding quantities entering 
in eq. (2.10) for the velocity are defined analogously. 

Below, we will investigate the consequences of the linear model, eq. (2.9,2.10), on the 
resulting approximations for the density (in both real and redshift space), as well as for the 
momentum and kinetic energy densities. All of those quantities are functionals of x(q,rj) 
and v(q,rj), and therefore we denote them as Q[x, v] as we did in the Introduction. We 
already wrote down <5[a;,i;] and 5 s [a;,u] in Section 2.1. There, we can also read off the LoS 
momentum density (divided by the CDM particle mass) given by T,\ , and the LoS kinetic 

energy density (divided by the CDM particle mass) given by T,\ /2. 

To build an approximation for Q, one has to calculate the displacement fields s z and 
S 2LPT and their time derivatives for a given realization, and scale them with the growth factor 
functions to the redshift of interest. This can be done with any 2LPT code for calculating 
the initial conditions for N-body simulations. One then multiplies (in Lagrangian Fourier 
space) the Lagrangian displacement/velocity fields with the displacement/velocity transfer 
functions to obtain s+ and v+. The transfer functions themselves can be obtained from a 
small set of realizations from eq. (2.11) and the analogous expression for the R v, s. 

The next step is to include the effects of the mode coupling terms. One possible way 
of dealing with that problem is to make a model for smc and vmc, e.g. by assuming they 
are Gaussian for simplicity. Another route (which works better as we will show later) for 
including their effects is to first calculate the approximation Q+(k, rj) = Q[q + s+, v+](k, rj) 
(where k can be a wavevector corresponding to x, x s or q depending on the quantity Q) and 
then write the following linear model for Q: 

Q{k, n) = R Q (k, ri)Q+(k, r{) + Q M c(k, rj) , with (2.12) 

(Q(M)Q*(M)> 



RQ(k,v) 



<Q*(M)Q*(M)> 



where again the transfer function Rq is defined in the least-squares sense, and Que is given 
by (Q+(k,ri)Q* MC (k,r])) = 0. In this paper, we follow this second route as it requires no 
assumptions on the statistics of smc and vmc- Moreover, as we will show, Q+ and Q are 
very well correlated, implying that Q+ is indeed a good approximation for Q. 

To summarize, the model we pursue in this paper is given by the equations of this 
section: eq. (2.9, 2.10, 2.12). Therefore, estimators (denoted by a hat) of fully non-linear 
quantities, Q, can be built for each realization according to the following scheme: 

{slpt, s LPT ) -> («•,«*) ->• Q* ->• Q , (2.13) 

once one has calculated the corresponding transfer functions. The same applies if one wants 
to construct unbiased estimators of statistics of Q such as its power spectrum. In that case, 
the mode coupling term, Qmc, can be dealt with in the same way as done in TZ. Note that 
since eq. (2.12) is a linear model for Q, it should be possible to invert the arrows in the 
scheme above, thus obtaining a self consistent reconstruction technique for the BAO peak. 
However, we postpone such a study to future work. 



3 Numerical results 

In this section we present numerical results based on the model given by eq. (2.9, 2.10, 2.12). 
The results that follow are extracted from 26 N-body simulations, made using GADGET- 
2 [14], with 256 3 particles each, in a box with side Lbox = 500Mpc//i. The cosmological 
parameters we used are (in the standard notation) as follows 

(fib, Matter, «a, K n„ <t 8 ) = (0.046, 0.28, 0.72, 0.70, 0.96, 0.82) . (3.1) 

The initial conditions were set at a redshift of 49, using the 2LPT code provided by [15]. 

Following TZ, we quantify the quality of the model given by eq. (2.9, 2.10, 2.12) by 
plotting the cross-correlation coefficient, pn, between a fully non-linear quantity Q and its 
approximation, Q+: 

<q(M)q;(M)> 

PQ V<IG*(M)I 2 XIQ(M)I 2 > ' 

where in this work Q stands for any of the following: s, v, 5, 6 S , k\\TS , k?,T,, , and the 
wavevector k above corresponds to q, q, x, x s , x, x, respectively. 

The plots that follow are done for two redshifts: z = and z = 1. Unless stated 
otherwise, for each quantity Q[a5,u] we show Rq and pq (eq's. (2.12, 3.2)) for four different 
approximations, Q+. Two of them correspond to the ZA and 2LPT results: 

Qir = Qz = Q [q + s z ,v z ] and Q+ = Q 2 lpt = Q[q + s z + s 2 lpt, v z + v 2 lpt] , (3.3) 

which also correspond to the model in TZ, where the particle positions and velocities were 
not corrected with transfer functions. The next two Q*'s are given by 8 : 

Q-k = Q[q + Rz*s z ,R v z *v z ] and (3.4) 

Q* = Q W + R z * $z + R2LPT * S 2 LPT, R" z * V z + R2LPT * V 2 LPt] , 

which correspond to the first and second order of the model of this paper, eq. (2.9, 2.10, 
2.12). The quality of any of the Q-^-'s at a given scale is measured by how close po(k, rj) is 
to one. 

It is important to stress that by analyzing all the approximations above, we do not 
aim to search for the best of them. It is by construction that the second approximation in 
eq. (3.4) must be the best. If any of the other approximations happen to be better for some of 
the quantities we analyze, that must be entirely for the wrong reasons. We will see examples 
of this as we proceed. Therefore, the comparison presented in the rest of this section should 
ultimately be treated in two ways: 1) As providing a better physical intuition of what effects 
are captured by which order in the LPT expansion (with or without transfer functions); 2) 
As giving us some idea of what improvements one may expect if one goes to higher order 
LPT. 



8 A natural extension to our models is calculating the 2LPT terms using an optimally smoothed 6l- How- 
ever, after applying such a smoothing, we did not find an appreciable improvement to the resulting non-linear 
5, which is why we stick to the simpler models above. 
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3.1 Displacement and velocity transfer functions 

In the top row of Figure 3 we show the velocity and displacement transfer functions 9 (R z , R v z , 
R-2LPT, ^2LPt)'i while the bottom two rows show the cross-correlation functions between s 
and its approximation s+, as well as between v and v+, for the four different approximations 
listed in eq. (3.3, 3.4), in which one should substitute Q with either s or v. 

Since at large scales the single-stream assumption behind LPT is not broken and the 
higher orders become negligible, all transfer functions and cross-correlations should go to 1 
for k/k^L — > 0, which they indeed do (Fig. 3). As can be seen from the red curves in the 
plots at the bottom of Fig. 3, at low k the cross-correlation coefficients between the ZA and 
the NL displacements and velocities approach 1 as 5 2 (< k) (defined below), which is the 
density variance integrated up to k. This is not surprising since the ZA in Lagrangian space 
misses second-order terms of order the large-scale density variance (see TZ). 

Including the 2LPT contributions clearly makes the low-A; power-law dependence much 
steeper, since 2LPT accounts for the terms in the transfer functions of order 5 2 {< k), and 
misses only higher-order contributions. This result will appear again and again considering 
other quantities, such as the CDM density, as well as the momentum density (the kinetic 
energy density is a special case, which we will consider in detail below). 

The different contributions to the displacement and velocity fields for our best model 
(second line of eq. 3.4) are shown in Fig. 4. The MC vector fields are split (in Lagrangian 
space) into irrotational and solenoidal parts. The solenoidal power is smaller than the irrota- 
tional power both for the velocity and the displacements, since virialization is not complete 
at these scales. The mode-coupling power for the velocity becomes comparable to the ap- 
proximated velocity, v+, at the non-linear scale (fc/vx = 0.25/i/Mpc), while this is not the 
case for the displacement vector field. This can easily be explained by looking at Fig. 3 (see 
also Section 3.1.1 below): At high k > kjsiL the LPT and NL quantities decorrelate since LPT 
breaks down after shell-crossing. The velocity cross-correlation and transfer functions decay 
faster than the corresponding quantities for the displacements, resulting in correspondingly 
larger velocity MC power. This result will be addressed next by constructing a simple model 
for the time dependence of R v z . 



9 Here we provide analytical fits (for the cosmology given in (3.f )) for the two displacement transfer functions 
and the density transfer function for our best density approximation (second line of eq. (3.4)): 

k - 2 



R z (k,rj) = exp -0.085 

«°™ m = exp (°- 64 (dy - l? yy J + °- 623 (dy J - °-° 78 y^ 



k N L(v) , 

2 / 7 \ 3 



Rs(k,rj) = exp(0.58d) , with d = S 2 ( < r^ k NL(r)o), no ) , (3.5) 

V k NL (ri) J 

with i]o evaluated at z = 0, and kNL and 5 2 can be read off from eq. (3.9). These equations were checked for 
2 = and z — 1 and give a respective accuracy of within (2%,4%,1%) at z = and (4%,4%,1%) at z = 1 
in the range 0.02/i/Mpc< k < 0.5/i/Mpc. The velocity transfer functions and the transfer function for the 
density in redshift space do not have a simple time dependence as above, and so we do not include fits for 
them. 
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Figure 3: Cross-correlation and transfer functions for the CDM displacements and velocities. 
In each plot we show the denoted quantities for the four different approximations given in 
eq. (3.3, 3.4). The red curves show the corresponding cross-correlation function divided by 
<5 2 (< k,r]), to highlight the low-k behavior. The non-linear scale (defined here as the linear 
power per logarithmic A:-bin to be 1) is at k^L = 0.25/i/Mpc for z = (denoted with vertical 
line) and k^L = 0.74/i/Mpc for z = 1. 
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Figure 4: In this figure we show the different contributions to the power spectrum of the 
displacements (left panel) and velocities (right panel). The total quantities are denoted with 
NL, and the mode-coupling vector fields are split (in Lagrangian space) into irrotational 
and solenoidal parts. The plot is for z = 0; and here k is the wave-vector corresponding to 
Lagrangian space. 

3.1.1 Why are the displacements better reconstructed than the velocities? 

In this section we offer a toy model relating the velocity and displacement transfer functions 
in the ZA. It should help us understand better the decay of the velocity transfer function, 
relative to the higher- k decay of the displacement transfer function, as observed in Fig. 3. 




0.02 



0.05 



o.io 

k[h/Mpc] 
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0.50 



Figure 5: Toy model for the ZA displacement and velocity transfer functions evaluated at 
z = 0. 



Starting from 



R z — d, 



{s z (t]) -s(r])) 
(\sM\ 2 ) 



dr, 



1 (s z {r] ) ■ s(r])) 

D(v) <M%)I 2 > 



(3.6) 
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and using s z (rj) = Ds z (rjo), we can write: 

R z = - (R v z - R z ) , with R z {k, n = 0) = 1 . (3.7) 

Since k^L is the relevant scale of the problem, we can write the time dependence for R^.(k, rj) 
as follows (see eq. (3.5) as well): 

Rlik^^R^k^M,^ , (3.8) 

where fc/vx(?7) is the solution to 10 

rkNLiv) 
S 2 (< k NL , rj) = 4tt / dkk 2 P L (k, r, )D 2 { v ) = 1 , (3.9) 

Jo 

and Pl is the linear power spectrum, defined as (<Jz,(fe, r))8L,(k, T))) — $D(k + k)Pi,(k,r]). 

As a check of the proposed time dependence of R v z , in Figure 5 we plot R v z for z = as 
obtained in two ways: directly from the N-body simulations (line denoted "vel. (ZA)"); as 
well as from (3.8) evaluated with rjo corresponding to z = 1, and rj to z = (line denoted 
"model vel. (ZA)"). We can see that there is a very good agreement between the two curves, 
which indicates that we capture most of the time- and scale-dependence of R v z with (3.8). 

Using (3.8), we can then integrate (3.7) to z = starting from the initial conditions at 
rj = 0, eq. (3.7). Thus, we obtain an R z at z = 0, which we plot in Figure 5 (line denoted 
"model displ. (ZA)"). This should be compared with the line denoted "displ. (ZA)", showing 
the true R z . Again we see a very good match between the two. So, we find that the time 
integral relating the velocity and the displacement moves the decay in the transfer function 
to higher k. This can be expected, since displacements are more insensitive to late time 
non-linearities than the velocity. To see this, one simply has to consider a virialized halo, 
where the velocities (and accelerations) can grow arbitrarily large, while the NL contribution 
to the displacement is bound to the size of the halo and is typically much smaller than the 
large-scale (linear) displacements. 11 



10 Elsewhere we quote fcjvi as the scale where the linear power per logarithmic fc-bin becomes 1. However, 
we find that (3.9) gives better agreement in eq. (3.8). 

The reader may be interested in the transfer function for the accelerations. So, let us look into it in some 
detail. Above we obtained the time dependence of R z to a good approximation by integrating (3.7). One can 
then easily calculate the transfer function, R s ii, relating the accelerations s" and s" . Using the properties 
of the Zel'dovich displacement as we did above, we can write R s » — {DR Z )" / D" . The result is the curve 
denoted "model s" (ZA)" in Figure 5. One then may be tempted to use the equation of motion (2.7) and 
write R s " ~ Rl, where Rl is the transfer function relating the linear and NL overdensities. This relation 
comes about because naively the proportionality (k- and time-dependent) constant between s" and S cancels 
out when calculating those transfer functions. 

However, this guess is not correct, since the k in R s n(k,ij) corresponds to Lagrangian Fourier space, while 
the k in RL{k,rf) is in Eulerian Fourier space. Thus, the large-scale dependence of R s ij(k,rj) can be only 
on the overdensity. That can be seen from eq. (2.7), where adding an overall constant vector to s leaves 
the equation invariant. Therefore, large-scale bulk motions are irrelevant for R s n . However, the large-scale 
dependence of RL(k,ri) (shown with the curve "density (linear)" in Fig. 5) depends on the large-scale flows 
as well, and thus on the large-scale velocity dispersion (see TZ). Therefore, it is not surprising that R s ij and 
Rl are qualitatively different - the former exhibiting a steeper decay than the latter. 
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Figure 6: Slices in real space, 100Mpc//i on the side (z = 0). 

3.1.2 Accuracy of the reconstructed particle positions. 

Let us now analyze the accuracy of the predicted positions and velocities, relative to the 
ones obtained by N-body simulations. To illustrate the effects of including the displacement 
and velocity transfer functions, in Figure 6 we plot slices, 100Mpc//i on the side, 7.8Mpc//i 
thick, through one of our N-body simulations in real space. The red dots denote the "true" 
particle positions as calculated from the N-body code. Each red dot is connected with a 
black segment to the approximated particle position (x = q + s+), obtained according to 
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one of the four approximations given in eq. (3.3 and 3.4), corresponding to the four different 
panels in each figure. 

One can confirm that even the ZA, the simplest approximation we use, reproduces the 
qualitative features of structure formation. We can also see that 2LPT over-corrects the 
particle positions in the voids (i.e. most black segments flip from one side to the other of the 
red dots when going from the upper left to the upper right panel of each figure) , while in the 
overdense regions we see that halos become more diffuse in 2LPT. After the displacement 
transfer functions are included (lower two panels), we can see that the voids are modeled 
almost perfectly at the resolution of the plots (as anticipated in the Introduction); and halos 
and filaments are modeled to a better accuracy than using pure LPT. 

To further quantify how well we reproduce the particle positions in real space, in the 
upper panels of Fig. 7 we plot the probability density for finding a volume element in Eulerian 
space quantified by (log(l + 5+) ,log(ArcM)) at z = 0, where both quantities are calculated 
after Gaussian smoothing on a scale of 1.3Mpc//i. Here ArcM is the error we make in 
predicting the center of mass (CM) position of that volume element by using the displacement 
approximations (given by eq. (3.3, 3.4)). In the middle row of Fig. 7, we show the same but 
for the probability density in the (log(l + 5^) ,log(AvcM/'H)) plane; while in the bottom row 
we show the probability density in the (log(l + 5^ r ),log(Av Tms /'H)) plane. This figure then is 
the 2D analog of Fig. 2, showing explicitly the density dependence. 

One can see that the 2LPT approximation including the displacement transfer function 
is the best at predicting the CM position both in underdense and overdense regions, the 
former being markedly better reproduced (especially at the 95% level) in agreement with our 
qualitative assessment from Figure 6. The high-density region of this plot can be used as a 
guide for the accuracy of our displacement approximations to reproduce halo positions both 
in real and in redshift space. Thus, we find that for our choice of smoothing, the position 
of most (looking at the 2-sigma, or 95% contour) overdense regions are found to better 
than ~3Mpc//i. Concentrating on the 1-sigma contours in the upper middle panel, where we 
compare the 2LPT results with and without the transfer functions, we find an improvement in 
the approximated CM positions which is about a factor of 2 better (corresponding < 1 Mpc//i) 
in real space. The FoG effects may introduce problems in redshift space as their effects are not 
captured as well: Looking at the bottom row, we can see that the errors in the rms velocity 
can be as large as lOMpc/h (when mapped to real space) in high density regions, which 
is comparable to the Doppler displacement from the velocity dispersion in galaxy clusters. 
However, the 1-sigma contour is again within a couple of Mpc, or a factor of a few smaller. 

3.2 Density in real space 

Next, in Fig. 8 we plot the cross-correlation and transfer functions for the density in real 
space for the four approximations of 5 given in eq. (3.3, 3.4) (with Q replaced by 5, which is 
given by eq. (2.2)). For comparison, we also include the transfer function Rl relating 5l with 
5, as well as their cross-correlation, pi- As argued in TZ, 5l is proportional to the initial 
overdensity, and therefore Rl and pi involve statistics of quantities at unequal times. Thus, 
for scales smaller than the rms particle displacements (k ~ 0.12/i/Mpc), the cross-correlation 
PL is destroyed (see TZ for detailed discussion). This is not the case for the 5+'s based on 
LPT, where the only relevant scale is Jznl (see TZ). 

Concentrating on the rest of the transfer and cross-correlation functions, we can notice 
several interesting results. We find that including the displacement transfer functions (the 
velocity is irrelevant for 5) as per eq. (3.4) has no effect on the large scales. That result is 
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Figure 7: Probability density in the (log(l + <5^-),log(ArcM)) plane (upper row), in the 
(log(l + <5*),log(A'UcM/'W)) plane (middle row), and in the (log(l + ^),log(Ai; rms /H)) plane 
(bottom row). Here, 5+ is the density in real space (obtained after Gaussian smoothing on 
a scale of 1.3Mpc//i) as predicted by the approximation denoted in the plot. The residual 
displacement (not captured by our approximations) of the center of mass (CM) of a given 
cell in Eulerian space is given by ArcM- The residual CM and rms velocities of that cell are 
given by Avqm and Av rms , respectively. For each approximation we show the 1- and 2-sigma 
contours with solid and dashed lines, respectively. The plot is for z = 0. For that redshift 
the rms linear displacement is « 15 Mpc/h, or about an order of magnitude larger than the 
residuals depicted in the plot. 
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Figure 8: The same as in Fig. 3 but for the density in real space. 

expected, because before shell-crossing LPT is exact and the displacement transfer functions 
go to 1 as discussed before. Another thing to note is that including both R z and R2LPT 
makes the low-A; R$ (top row of Fig. 8) flatter, lying closer to 1. 

At small scales, we improve significantly the cross-correlation between 6+ and 5. Thus, 
a level of p$ = 0.95, is reached at k ~ 0.37/i/Mpc for the 2LPT approximation (the best ap- 
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proximation used in TZ), while including the displacement transfer functions to second order 
pushes that scale to k ~ 0.45/i/Mpc. Note that without the displacement transfer functions, 
the ZA and 2LPT cross-correlation curves intersect in the NL regime, while this is no longer 
the case after including the displacement transfer functions. This result can be understood 
by the fact that R z and R2LPT capture some of the higher-order LPT contributions (missed 
by our second order truncation), as well as some of the stream-crossing effects missed by 
LPT altogether. 

In TZ we showed how one can build an approximation for 5 using LPT to reduce the 
sample variance of the matter power spectrum. The amount of the sample variance reduction 
is approximately given by the ratio of Rj(\5+\ 2 ) to the mode coupling power, (\5mc\ 2 )- This 
boils down to Pg/(1 — pf) which for p$ ~ 1 is simply ~ 1/(1 — pf). This sample variance 
reduction corresponds to a speed-up of the scanning of the cosmological parameter space one 
can gain. Thus, one should pay close attention to the bottom row of Fig. 8, showing 1 — p 2 . 
We can see that by using our best 5+ to build an estimator for the matter power spectrum 
as done in TZ, we can gain a speed-up of a factor of 10 5 at scales (k ~ 0.025/i/Mpc) ten 
times larger than the non-linear scale; while at the non-linear scale the expected speed-up is 
about a factor of ten to a hundred. Therefore, we see that this speed-up is significant well 
into the NL regime, especially if one uses the transfer-function-corrected displacements. 

3.3 LoS momentum and kinetic energy densities 

Let us now focus again on the expansion of S s given by eq. (2.4). We already have a plot 
(Fig. 8) of the transfer and cross-correlation functions for the zero order (in the velocity) S s , 
which is simply 5 in real space. The properties of the next two orders are plotted in Figures 9 
and 10. For these plots, we make the flat sky approximation for simplicity. 

As one can see by comparing Figures 8 and 9, the density and the momentum ap- 
proximations share the same features. 12 That is, apart from the fact that all momentum 
transfer and cross-correlation functions deviate from 1 at larger scales - a direct consequence 
of the fact that the velocity transfer/cross-correlation function decays at larger scales than 
the corresponding displacement functions (see Section 3.1). 

From Fig. 10, we can see that qualitatively new features appear for 5 S at second order 
(in the velocity), which corresponds to the second LoS derivative of the LoS kinetic energy, 
k?,T,\ . Most strikingly, none of the p's for this quantity saturate to 1 at large scales. 

Moreover, the R and R" corrected approximation for k?,T,\ is well below the uncorrected 
2LPT result (at least at large scales for z = 0). The reason behind both of these results 
is that the kinetic energy even at large scales is sensitive to the high-/c velocity dispersion. 
This is because short-scale random motions (inside halos for example) do not cancel for the 
kinetic energy, as they do for the momentum where only CM motions matter at large scales. 
Since the uncorrected 2LPT result has a larger (albeit artificially larger) velocity dispersion 
than the result corrected with R v < 1, the kinetic energy coming from the short-scales is 
boosted, thus compensating for the missing short-scale power due to the missing vmc- 

To lowest order in Eulerian perturbation theory, one can check that the high-A; contri- 



12 The corresponding plots for the LoS momentum density approximations are very close to the ones for the 
LoS derivative of the LoS momentum density, Fig. 9, which is why we do not include the former. 
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Figure 9: The same as in Fig. 3 but for fciiTI| which is proportional to the LoS derivative 
of the LoS momentum density: dpu/dxu. 



bution to k?,T,\ goes as 



knT^ ~ fen / dkPv(k)k 2 (high k contribution only) , 



(3.10) 
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Figure 10: The same as in Fig. 3 (dropping the bottom row) but for k?,Tn which is 
proportional to the LoS second derivative of the LoS kinetic energy density: cPEj.\\/dx?, 

where P v is the velocity power spectrum. Therefore, for a given velocity approximation, 

1 — p 2 2 (2) can be estimated by the fractional difference between the quantity in eq. (3.10) 

fc ll T ll 
evaluated with the fully non-linear P v , and the quantity in eq. (3.10) evaluated with the P v 

obtained from the said velocity approximation. We find that at small k this fractional differ- 
ence is indeed on the order of several percent for all approximations we consider, consistent 
with the bottom row of Fig. 10. 

Since k?,T,\ contributes to S s , we can expect that p$ s will not be as well reconstructed 
as one might expect from the ps plots (Fig. 8). We will see that that is indeed the case in 
the next Section. 

3.4 Density in redshift space 

Now let us turn to Fig. 11 where we plot the same quantities as before but for the den- 
sity monopole in redshift space for the four approximations for 6 S given in eq. (3.3, 3.4) 
(with Q replaced by S s , which is given by eq. (2.3)). For these plots, we make the flat sky 
approximation for simplicity. 
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Figure 11: The same as in Fig. 3 but for the density monopole in redshift space. 

Comparing Fig. 8 and Fig. 11, we find that including the displacement and velocity 
transfer functions does not improve ps s significantly at z = 0, while at z = 1, the presence of 
the transfer functions makes p$ 3 even worse at high k. As we already showed in detail in the 
previous section, this effect is due entirely to the wrong reasons: Artificially setting R v = 1 
captures by accident some of the effects of the mode-coupling velocities, which are missed by 
our best approximation (second line of eq. (3.4)), for which R v < 1. We should stress that it 
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is completely unclear whether this result will persist for halos or other tracers of the density 
field. 

Nevertheless, let us try to understand this result in more detail. We should remember 
that FoG effects move power from the NL to the MNL regime. By including the R v, s, which 
decay at high k, we reduce the high-/c velocity power, which makes the FoG effect smaller. 
Thus, because of the non-linear mapping, eq. (2.2), between displacements and velocities on 
the one side and 5 S on the other, we find that using R v < 1 at high k, affects the MNL 5 S .+ . 
Those effects cannot be captured by the linear (isotropic) rescaling, R$ s , and are instead 
dumped into S St MC- 

To test that this effect on ps a is indeed due only to the R"'s being smaller than one, 
we constructed a d St + by including only the R v, s (i.e. setting R z = R2LPT = 1) and then 
including only the displacement transfer functions (i.e. setting R" = R^lpt = -0- We ^ n< ^ 
that the former choice results in degradation of ps s , while the latter choice improves p$ s for 
both redshifts compared with the pure 2LPT result. Thus, we find that faking the high-/c 
power in vmc by keeping R v = 1 improves ps s - However, such a route is not consistent, 
since the mode-coupling and LPT velocities must have a vanishing 2-pt function. Instead, 
if one is determined to boost ps a closer to 1, one may be better off including a Gaussian 
vmc-, uncorrelated with s z / 2 lpt- However, this would result in adding a velocity dispersion 
both in halos and in voids. Nevertheless, we tested that approach 13 , and indeed it results in 
a reduced ps 3 , leading us to conclude that the higher order correlations between vmc and 
s z /2lpt are important. Indeed looking at the bottom row of Fig. 7, we can see that the rms 
residual velocity is indeed correlated with the overdensity. However, further analysis along 
those lines will be warranted only if one calculates these transfer function for CDM halos 
(not the CDM matter density) or other biased tracers. 

4 Summary 

In this paper, we consider an improvement to the formalism proposed in TZ for studying 
the mildly nondinear regime of structure formation. This improvement is obtained by con- 
structing good approximations of the CDM particle trajectories. This is done by relating 
the displacement and velocity fields obtained in LPT (in Lagrangian space) to the non-linear 
ones by using time- and scale-dependent transfer functions. We work up to second order in 
LPT. Compared to the pure 2LPT results, after including the transfer functions we find an 
improvement in the approximated center of mass position of a given cell (in Eulerian space) 
which is about a factor of 2 better (< 1 Mpc/h). 

These improvements in approximating the particle positions result in an improved ap- 
proximation, 5+, for the fully non- linear density, 5. The cross-correlation coefficient between 
5+ and 5 reaches a level of p$ = 0.95 at k ~ 0.37/i/Mpc for the 2LPT approximation (the best 
approximation used in TZ), while for our current best 5+, which includes the displacement 
transfer functions to second order, the corresponding scale is k ~ 0.45/i/Mpc (for z = 0), 
almost twice smaller than the non-linear scale. For comparison, the corresponding scale for 
the cross-correlation between the density from linear theory and 5 is k ~ 0.075/i/Mpc, a fac- 
tor of 6 larger, than our best density approximation. Thus, our results allow for a speed-up 
of an order of magnitude or more in the scanning of the cosmological parameter space with 
N-body simulations for the scales relevant for the baryon acoustic oscillations. 



13 We added vmc by adding both its solenoidal and irrotational components shown in Fig. 4. Doing the 
same for the displacements only further worsens the results. 
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We also investigated the effect on using our trajectory approximations for calculating 
the density in redshift space. We find that compared to S s obtained in pure 2LPT, the 
resulting 8 S approximation has only slightly improved cross-correlation with the NL 5 S for 
z = 0. That cross-correlation is even slightly worse than the 2LPT result for z = 1. However, 
we argue that this slightly better cross-correlation for 2LPT is due to an inconsistency in the 
2LPT approximation, and therefore one should still use our best approximation. Moreover, 
it is unclear whether this result will hold for halos or other biased tracers - an analysis which 
we postpone to future work. 

As an application of the results presented here, in a companion paper [16], we utilize the 
models of the non-linear density field constructed in this paper, to develop a quasi-optimal 
reconstruction scheme of the BAO peak. 
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